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Abstract 

The following questions are often encountered in system and control 
theory. Given an algebraic model of a physical process, which variables can 
be, in theory, deduced from the input-output behavior of an experiment? 
How many of the remaining variables should we assume to be known in 
order to determine all the others? These questions are parts of the local 
algebraic observability problem which is concerned with the existence of 
a non trivial Lie subalgebra of the symmetries of the model letting the 
inputs and the outputs invariant. 

We present a probabilistic seminumerical algorithm that proposes a 
solution to this problem in polynomial time. A bound for the necessary 
number of arithmetic operations on the rational field is presented. This 
bound is polynomial in the complexity of evaluation of the model and 
in the number of variables. Furthermore, we show that the size of the 
integers involved in the computations is polynomial in the number of 
variables and in the degree of the differential system. 

Last, we estimate the probability of success of our algorithm and we 
present some benchmarks from our Maple implementation. 

Keywords: Local algebraic observability, local algebraic identifiability, 
seminumerical algorithm. 

Mathematics Subject Classification (2000): 93B07, 93B40, 93A30; 12H05. 



1 Introduction, Notations and Main Result 

Local algebraic observability is a structural property of a model and one of the 
key-concepts in control theory. Its earliest definition goes back to the work of 
R.E. Kalman for the linear case (see ^]J) and a large literature is devoted to 

"This paper is available at B8I. All comments are welcome. 
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figure 1: Model for circadian oscillations in the Drosophila period protein jf?] 
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this subject (see £l6|, fll], [Tl| and the references therein). We base our work 
on the definition given by S. Diop & M. Flicss in jllj of the observability for 
the class of algebraic systems. 

As in the example of figure [l], such a system is usually described by means of 

• a vector field, which describes the evolution of state variables in function 
of inputs and of parameters; 

• some outputs which are algebraic functions of these variables. 

The definition of observability given in jll] relics on the theory of differential 
algebra founded by J.F. Ritt Q and is based on the existence of algebraic 
relations between the state variables and the successive derivatives of the inputs 
and the outputs. 

These relations can be considered as an obstruction to the existence of in- 
finitely many trajectories of the state variables which are solutions of the vector 
field and fit the same specified input-output behavior. If there are only finitely 
many such trajectories, the state variables are said to be locally observable. 

In order to illustrate this notion, let us consider the local structural identi- 
fiability problem which is a particular case of the observability problem. The 
question is to decide if some unknown parameters of a model are observable con- 
sidering their parameters as a special kind of state variables satisfying O = 
(see Q Q §). If they are not observable, then infinitely many values 

of these parameters can fit the same observed data. Hence, if these parameters 
have a physical significance, it may be necessary to change the experimental 
protocol when possible. On the other hand, if the parameters are identifiable, 
various numerical approximation methods can be used for their estimation (see 
| jj9| and the references therein). 

We consider the local algebraic observability problem under the computer al- 
gebra standpoint. The previous studies that enable to test observability mainly 
rely on characteristic set or standard bases computation |33[ ^9], ^, |l^] and 
their complexity is, at least, exponential in the number of variables and of pa- 
rameters (see [ful 5||). Some other techniques, as the local state-space isomor- 
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phism approach E3] or the conversion between characteristic set w.r.t. different 
ranking ||, can also be used. The complexities of these methods are not known. 

We present a probabilistic polynomial-time algorithm which computes the 
set of observable variables of a model and gives the number of non observable 
variables which should be assumed to be known in order to obtain an observable 
system. A Maple implementation is available at |38| . 

Example: Let us consider the use of our algorithm with a model for circadian 
oscillations in the Drosophila period protein []l7|| . This model is presented in 
figure [j]; there are seventeen parameters and no input in it. After one minute 
of computation, our Maple implementation gives the following results: 

• the variable M and the parameters {v s , v m , K m , k s } are not observable. 
All the other parameters and variables are observable; 

• if the non observable variable or only one of the non observable parameters 
are specified, all the variables and parameters of the resulting system are 
observable. 

Our algorithm certifies that a variable is observable and the answer for a 
non observable one is probabilistic with high probability of success. These results 
allow us to focus our attention on just four of the seventeen original parameters. 
Thus, the search of an infinitesimal transformation which leaves the output y 
and the vector field invariant is simplified and we find a group of symmetries 
generated by {M,v s ,v m , K m , k s } — > {AM, Xv s , Xv m , XK m , k s /X}. Hence, there 
is an infinite number of possible values for non observable parameters which fit 
the same specified output y: this system is certainly unidentifiable. 

1.1 Notations and Main Result 

Hereafter, we consider a state-space representation with time invariant param- 
eters defined by an algebraic system of the following kind: 



Big letters stand for vector-valued objects and we suppose that there are: 

• I parameters :— (6\, . . . , Oi) 

• n state variables X := (x\, . . . , x n ); 

• r input variables U := (ui, . . . , u r ); 

• m outputs variables Y := (yi, . . . , y m ) with m < n. 

The letter X stands for the derivatives of the state variables (±1, . . . , x n ) and F 
(resp. G) represents n (resp. m) rational fractions in Q(X, 9, U) which are 



E 




0, 

F(x,e,u), 

G(X,Q,U). 



(1.1) 
(1.2) 
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denoted by , f n ) (resp. (ffi, ■ ■ ■ ,3m))- The letter d (resp. h) represents 

a bound on the degree (resp. size of the coefficients) of the numerators and 
denominators of the /j's and <j>j's. 

Hereafter, we use a common encoding where the expression e := x 5 is repre- 
sented as a sequence of instructions: t% :— x, t 2 := ii 2 , £3 := t 2 2 , e :— t$t\. 

Hence, the system E is represented by a straight-line program without divi- 
sion which computes its numerators and denominators and requires L arithmetic 



operations (see Section |3.4| and § 4 in |]9j) 



The following theorem is the main result of this paper. 



Theorem 1 Let £ be a differential system as described in Section 1.1. There 
exists a probabilistic algorithm which determines the set of observable variables 
of £ and gives the number of non observable variables which should be assumed 
to be known in order to obtain an observable system. 

The arithmetic complexity of this algorithm is bounded by 

O ]M{v) (N{n + £) + (n + m)L \ + (n + £ + l)N(n + I) 



n + £l 

with M(y) {resp. N{y)) the cost of power series multiplication at order v + 1 
(resp. v x v matrix multiplication) where v < n + £. 

Let [i be an arbitrary positive integer, D be 4(n + £) 2 (n + m)d and 

D' := (2 ln(n + £ + r + 1) + In uD) D + 4(n + £) 2 ((n + m)h + In 2nD) . 

If the computations are done modulo a prime number p > 2D' '/i then the proba- 
bility of a correct answer is at least (1 — l/fi) 2 . 

For the model presented in figure [l], the significant terms of our complexity 
statement are L = 91, n = 5,£ = 17, m = 1, d = 6, h = 1, v = n + 1. The choice 
of fi = 3000 leads to a probability of success around .9993 and the computations 
are done modulo 10859887151. These computations take 10 seconds on a PC 
Pentium III (633 Mhz) provided by the UMS MEDICIS @. 



Outline of the paper: In the next section, we recall some basic defini- 
tions of differential algebra and the definition of algebraic observability used 
by S. Diop & M. Fliess in (T^|. Furthermore, we describe the relationship be- 
tween this framework and the approach of H. Pohjanpalo in |3^|. Then, we 
present an algebraic jacobian matrix which is derived from the theory of Kahler 
differentials and used in the local algebraic observability test. 

In the second part of this paper, we present some new results. In Section [| 
we show how to compute some specializations of this matrix using power series 
expansion of the output and we estimate the related arithmetic complexity. 
Then, we study the behavior of the integers involved in the computations and 
we precise the probabilistic aspect. In conclusion, we present some benchmarks. 
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2 Differential Algebra and Observability 



Differential algebra, founded by J.F. Ritt, is an appropriate framework for the 
definition of algebraic observability introduced by S. Diop & M. Fliess in pj| . 
For more details on differential algebra, we refer to [|4| and Hj; nevertheless, 
we recall briefly some necessary notions. 

2.1 Differential Algebraic setting 

Let us denote by k a base field of characteristic zero. The differential alge- 
bra k{U} is the fc-algebra of multivariate polynomials defined by the infinite 
set of indeterminates {£/'•-'' |Vj G N*} and equipped with a derivation 8 such 
that Su^ = u^ l+1 \ Its differential fraction field is denoted by k(U). 

Hypotheses: The inputs U and all their derivatives are assumed to be inde- 
pendent. Furthermore, we consider non singular solutions of S; thus, we assume 
that we work in an open set where the denominators present in E do not vanish. 
These hypotheses represent practically all the encountered systems. 



2.2 Local Algebraic Observability 

Following the interpretation due to M. Fliess of some algebraic control theory 
problems fll^l , we consider the differential field K. := k(U)(X, 9) equipped with 
the following formal Lie derivation: 



This derivation is associated with the vector field defined by the equations (1.1). 
Hereafter, we denote . . . , Cf n ) by CF and £ o ■ ■ ■ o £ by D . 

j times 

Hence, the outputs G(X, 6, U) are denoted by Y and = &G(X, 6, U). 



Definition 1 ([26, 11]) An element z in K, is locally algebraically observable 
with respect to inputs and outputs if it is algebraic over k(U,Y). Thus, the sys- 
tem E is locally observable if the field extension k(U, Y) > K. is purely algebraic. 

Let us illustrate this definition with the following example: 



£3 


= 0x\, 




= x 3 /x 2 


X\ 


= x 2 /xi 


V 


= Xi. 



By successive differentiations of the output, we obtain the following differential 
relations: 

y-xi, yy-x 2 , yy(y 2 + yy)-x 3 
(y 2 + yy) 2 + yy{3yy + yy {z) ) - Oy- 
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Thus, the parameter and the variables are observable according to Definition [l]. 
Furthermore, as these algebraic relations define a unique solution, the parameter 
and the variables are said to be globally algebraically observable [^6[ ^9], |9| . 

These relations depend generically of high order derivatives of the output 
and thus, they are not of a great practical interest for parameter estimation. 
As we focus our attention on local observability, we are going to avoid their 
computation. 

Convention: We wish to test observability for the parameters and/or state 
variables X . Thus, we present the algorithm in the most general case (observ- 
ability of parameters and state variables) and we do not describe the restriction 
to one case or the other. 

Definition |l| implies that local algebraic observability is related to the tran- 
scendence degree of the field extension k(U, Y) <^-> JC. So, this proper ty c an be 
tested by a rank computation using Kahler differentials (see Section 2.4). As 
noticed in |Tl} |, this approach leads to a condition which is the formal counter- 
part of the R. Hermann & A. Krener rank condition in the differential geometric 
point of view jl8| . 

Furthermore, the transcendence degree of the field extension k(U, Y) <—> JC is 
the number of non observable variables which should be assumed to be known 
in order to obtain an observable system. Thus, Theorem [l] is based on the study 
of this field extension. 



2.3 A Description of k(U, Y) ^ K 

Let us denote by $(A, Q,U, t) the formal power series with coefficients in JC 
such that <J>(A, 6, U, 0) := X and $ = 9, U), we have: 

<t>(X,Q,U,t) = X + ]T aF{X,Q,U) - v 

Furthermore, let us define the formal power series Y(X, 0, U, t) with coefficients 
in JC such that Y(X,@,U,t) := G($(X,Q,U,t),e,U,t): 

Y(x, e, u, t) = g(x, e,u)+J2 &g(x, e,u)^. (2) 

We recall that these expressions are vector- valued (Y = {y\, . . . ,y m )). 

In p]] ]) H. Pohjanpalo already considers the coefficients of the power se- 
ries Y(X, <d,U,t) in order to test identifiability. In pT| , the authors prove that 
a finite number of these coefficients are necessary to describe the field exten- 
sion k(U, Y) <^-> /C. But in these two papers the necessary order of derivation is 
not bounded. 

This can be done using the differential algebra point of view (see § 4 in |35| 
for a general statement). The following proposition summarizes these results in 
a field extension framework. 
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Proposition 1 The field k{U,Y) is isomorphic to k(U)(Y, . . . ,Y<- n+t+r >) and 
algebraic over k(U) (Y, . . . , Y~(™+0) . 

Proof: The transcendence degree of k(U) <— > K. is equal to n + t. Hence, the 
transcendence degree of k(U) °-> k(U,Y) is bounded by n + 1. It means that, 
for i = 1, . . . ,m, there is an algebraic relation qi(y%, ■ . . , yi" n+t> ) = and the 
derivative is a rational function of yi, . . . ,y^ n+ ^ with coefficients 

in k(U). This proves that fc(J7, V) is isomorphic to fc(C7) (Y, . . . , y(™+^>) . ■ 

If there is more than a single output, the necessary order of derivation can 
be smaller than n + 1 and it is denoted by v. This index of differentiation 
is a natural measure of the complexity of our algorithm (see Section 3.4) and 
generically v = (n + £)/m. Hereafter, we take v equal to n + £ as in Theorem 



In the above proof, following the hypotheses of Section 2T, we assumed that 
the independent input variables U and all their derivatives were in the base 
field. Furthermore, we showed that we just need the first n + £ derivatives of 
the output equations. In order to simplify the presentation in the next section, 
we assume that the base field is k := k(U, Y,..., U {n+e \ Y( n +Q) . 

We present now the properties of the module of Kahler differentials which 
are used to compute the transcendence degree of k k(X, 0) in practice. 



2.4 Rank Conditions 

If S » T is a field extension, we use the notation flx/s f° r the T- vector space 
which is the cokernel of the jacobian matrix d(£ l G) <i<v/d(X, 9) and dz stands 
for the image of z G T in this vector space (see § 16 in [Q for standard definition 
and [pil for construction in differential algebra). We recall the following result: 



Theorem 2 (§ 16 in ]12[ ) Let us consider S a field of characteristic zero 
and T a finitely generated field extension of S. If {xaIasA <ZT is a collec- 
tion of elements, then {Ax\\\^^ is a basis offlq^/s as a vector space over T iff 
the {xaIagA form a transcendence basis of T over S. 

Our algorithm is based on the following straightforward consequences of this 
theorem. 

Corollary 1 If '<p is the transcendence degree of the field extension k <— > k(X, 0) 
then we have the equality 

= (n + i) _ rank S(x , e) {d(CG) Q ^ u /d(X, 6)) . 

Furthermore, If the rank of the jacobian submatrix d(C 3 G)a<j< u /d(X\{xi},@) 
(resp. G)o<j<„ / d(X, ®\{6i}) ) is equal to n + i — <f>, then the transcendence 
degree of the field extension k > k(xi) (resp. k k{9i)) is equal to zero and 
the variable Xi ( resp. the parameter 9i ) is observable. 
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The computation of tf> is mainly based on the construction and the evaluations 
of a straight-line program which allows to compute the power series expansion 
of Y(X, Q,U,t). We present the necessary notions in the next section. 



2.5 Data Encoding and Complexity Model 

The above results can be expressed considering a polynomial / as an element 
of a vector space; hereafter, we consider an algebraic expression as a function. 

This classical point of view in numerical analysis is also used in computer 
algebra for complexity statements or practical algorithms (see |L6], [l5| p7| 
and the references therein). We refer to Chapter 4 of Q for more details about 
this model of computation. 

Definition 2 Let A := {ai, . . . , aj} be a finite set of variables. A straight-line 
program is a sequence of assignments bi <— b' Oj b" where o.j £ {+, — , x, -=-} and 
where {6', &"} C U}=i{^j} U AU k. Its complexity of evaluation is measured by 
its length L, which is the number of its arithmetic operations. Hereafter, we use 
the abbreviation SLP for straight-line program. 

As a SLP representing a rational expression / € k(a\, . . . , aj) is a program 
which computes the value of / from any values of the base field such that every 
division of the program is possible. Furthermore, it is possible to determine a 
SLP representing the gradient of /. The following constructive results allows us 
to handle these two aspects. 



Theorem 3 (W. Baur &; V. Strassen |l[ ) Let us consider a SLP computing 
the value of a rational expression f in a point of the base field and let us denote 
by Lf its complexity of evaluation. 

One can construct a SLP of length 5Lf which computes the value o/grad(/). 

Furthermore, one can construct a SLP of length 4L f which computes two poly- 
nomials /i and fi such that / = f\jfi- 

Following our presentation, one can construct formally all the expressions 



introduced in Sections 2^3 and 2A with its favourite computer algebra system. 

But, let us recall that, in order to compute the formal expressions C V G 
and the associated jacobian matrix, one has to differentiate v times the output 
equations (1.2). As explained in p2f , the arithmetic complexity of computing 
multiple partial derivatives is likely exponential in v. If the evaluation complex- 
ity of the output equations (1.2) is L, by Theorem ||, the computation of C V G 
requires at least (5m)"L arithmetic operations. 

Thus, this strategy cannot lead to a polynomial time algorithm. 

The rank computations defined in the previous section are also cumbersome 
because they are mainly performed on the field k(X, 0). Nevertheless, in order 
to determine (j> efficiently, the variables X, O and U can be specialized to some 
generic values in the jacobian matrix and so, its generic rank can be computed 



numerically with high probability of success (see Section 3.6) 
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Thus, the main problem is to avoid the formal computation of (DG)o<i<v 
In fact, our strategy is to specialize a linearized system derived form £ first and 
to recover the value of 4> just using numerical computations on a finite field. 

3 A Probabilistic Polynomial-Time Algorithm 

In Section [O], we present the linear variational system derived from £ which al- 
lows us to compute directly the jacobian matrix d(£ l G)o<i<u/d(X 7 0) with X, 
and U specialized on some given values. 

Then, we show how this matrix can be determined in polynomial time and 
we give an estimation of the arithmetic complexity of our algorithm. 



The purpose of the Sections 3.5 and 3.£ is to study the growth of the integers 



involved in the computations and to estimate the probability of success of our 
algorithm. 

3.1 Variational System Derived From S 



As shown in Section 2A, our goal is to compute the generic rank of the jacobian 
matrix d(£ l G)o<i<v / d(X , 6). Using relation (g), we conclude that: 

d(£ j G) Q < j < v _ d(cocffs(y(t))) _ fdG_d§_ 8G_ 9$ dG 

d(x, e) ~ d(x, e) ~ coeffs \dx dx ' ox do + oe 

The above equalities leads to the following relation: 

where W denote the following n x (n + 1) matrix represented by a SLP: 

vy(*XA,e, U): =(§r.f§A + ||(*,r,A,e, C ,). 

Hence, we have to determine the first v = n + 1 terms of the power series ex- 
pansion of 6, U, t), T(X, 8, U, t) := d^/dX and A(X, 6, U, t) := d<5>/dQ. 

Let us denote by P(X,X,Q,U) = 0, the numerators of the rational rela- 
tions X — F(X, Q, U) = and let us consider the following expressions: 

r p(x,x,q,u), (|i) 

vp I zz(x,e,u)t + §%(x,x,e,u)T, (I2) (4) 

{ §z(x,e,u)A+j£(x,x,e,u)A+%(x,x,e,u). §3) 

The power series <b(X,Q,U,i), T(X,Q,U,t) and A(X, Q, U,t) are solutions of 
the system of ordinary differential equations VP($, T, A, O, U) = with initial 
conditions T(X, 9, U, 0) := Id„ x „ and A(X, 6, U, 0) := 0„ x£ . 
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Commentary: We have already noticed that one can compute symbolically 
the expression of the formal jacobian matrix d(C l G)o<i<v / d(X , 0). The rank 
computations described in Corollary [l] are sufficient to conclude. 

Furthermore, if X, and U are specialized on some random values, these 
computations can be performed numerically with high probability of success. We 
summarize this possible strategy in the upper horizontal and the right vertical 
arrow of the following diagram: 



VP 



formal computation 



numerical computation on 



d(X,B) 



X 

e 

u 



\ >0<i 

d(X,0) 



■* x ez", 
-> eez ! , 

* u e (my 
{x ,e,u) 



As the symbolic computation of the jacobian matrix is cumbersome, we special- 
ize the parameters on some random integers and the inputs U on the power 
series U which are truncated at order n + £ + 1 with random integer coefficients. 
Then, we solve the associated system VP for some integer initial conditions Xq 
and we compute the specialization d(£ z G)o<i< 1J /d(X, <d)(X , 6) with VY. This 
approach is summarized by the left vertical and the lower horizontal arrow. We 
present an algorithm which relies on this standpoint and we give in Section |3 . 6| 
its probability of success. 



The hypothesis dP/dX ^ assumed in Section 2.1 ensures that the differen- 
tial system VP(<£>, T, A, 0, f/) =0 admits an unique formal solution |l^] which 
can be computed with the following Newton operator. 



3.2 A Quadratic Newton Operator 

The aim of this section is to present the Newton operator used in our algorithm. 
In |lj, ||, the authors show that such an operator is quadratic. We sketch its 
construction and neglect the technical details for the sake of simplicity. 

We recall that we work with vector- valued expressions. Thus, the expres- 
sion (Q.l) (resp. (4.2), (Q.3)) represents a n x 1 (resp. n x n, n x t) matrix. 

The Theorem 3 allows to construct, from a SLP of length L which encodes S, 
another SLP of length 0(N(n + £) + nL) which encodes the system VP. For 
some given series $,r and A, this SLP computes the following n x (1 +n + £) 
matrix: 

/ Pl ($,$,e,u) f!($,0,c7)r i#( $ '^) A +\ 

: + f£($,$,0,LOA + 

V vn (M,e,L0 |f (M.e.Lr) r fg(4>, <s>,e,u) J 
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Let us represent $(X,Q,U,t) (resp. A(X,G, U,t), T(X,G,U,t)) mod t 2 ' 
by $j (resp. Aj, Tj) and denote the correction term: 

(&(X,Q,U,t) -<S>3,T(X,G,U,t) - Tj,A(X,G,U,t) - Aj) mod t 23+1 by 

As usually, we construct our Newton operator from the Taylor series expan- 
sion of the function VP. This yields the following relations: 

3 yp . 8 VP 

vp($, r, a) (x, e, u, t) = vp(<s> 3 , r, , a,) + p E ^+ d (x, r, X) " = a 

The remaining terms are of order in t greater than 2 J+ . Thus, they are not 
necessary for the computation of Ej . 

Computational strategy: we consider $ as a variable in the first column 
of VP and as a constant in the others. Thus, we have the following relations: 

3 VP _ f 3P_ dP_ 3P_\ 8 VP _ f dP_ dP_ dP \ 
d(X,t,A) ~ \dX~ , dX~ , dX~)' d(X,T,A) ~ {dx'dx'dxj- 

Consequence of our computational strategy: The above hypothesis in- 
duces a shift between the order of correct coefficients of Aj, and <F,. In 
fact, Aj and Tj are correct modulo t 2J . Thus, we need to stop the follow- 
ing operator with j + 1 = lna(^ + (■ + 1) & n d to repeat one more time the last 
resolution at the same order. 

Newton operator: The above hypothesis leads to a Newton operator based on 
the resolution of the following system of linear ordinary differential equations: 

BP ~ ~ BP ~ ~ ~ ~ +i 

^ (*j , e, u) E j+1 + — fa, e, u)e j+1 + vp($ f , r,-, ^ ■., &, u) = o mod t 2 (5) 

From the initial conditions <&o G ? := Idnxn an d A := nx e, this system 
is solved iteratively for j + 1 = 1, . . . , ln 2 (n + £ + 1) using the recurrence rela- 
tions ($ J+1 ,r J+1 ,A J+1 ) = (Sj.r^Aj) +E j+1 . 

The resolution of the linear ordinary differential system (|5|) relies on the 
method of integrating factors. First, we consider the Homogeneous system 

BP ~ ~ BP ~ ~ +i 

^ (*,- , e, u) Wj + gx(i>i, *i, ©, u)w s = o mod 

where Wj denote a n X n unknown matrix which coefficients are series trun- 
cated at order 2 3 . The main trick is common in power series manipulation, 
we consider matrices with coefficients in a series ring as series with coefficients 
in a matrix ring. For example, we have A mod t 2 ' — Aq + Ait + • • ■ + A^jt 2,3 
where the A^s are matrices with coefficients in the rational field. 

Thus, the product, the exponential and, if Aq is invertible, the inverse of 
matrices with coefficients in a series ring can be computed at precision j with 
the classical Newton operator (see 4.7 in |2j| and § 5.2 in for more details). 
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figure 2: Local Algebraic Observability Test 



Input : X-F(X,Q,U), Y~G(X,<d,U) 
Output : Succeed, a boolean 

Preprocessing Construction of the SLP coding J^, VP, Re- 
initialization Choice of a prime number; 

U <— Random Power Series mod t n+e+1 ; 
Succeeds true; Order<— 1; Q<— Random Integers; 
A <— Onxe; r <— Id„xn; X <— Random Integers; 

while Order < n + £ + 1 do 

W <- HomogeneousResolution (f£($, O) W" + f^OX", $, 9) W = u)mod t' 

($, A, T) «- ($, A, r) + ConstantsVariation fw, VP ($, I\ A)) mod r° rdcr ; 

Increase Order; (Order <— 2 Order); 
end while 

JacobianMatrix <- Coeffs(VY ($, T, A) , t J j = 0, . . . , n + I) ; 

Test if n + l> Ran k( JacobianMatrix) 

then Succeed := false 
end if 



For example, if Aq is invertible and B 3 denotes the inverse of A at order t 2 ' , we 
have B j+ i = 2Bj - B 3 AB 3 . 

Furthermore, it is a basic fact from the theory of linear ordinary system that 
if AW + A'W = and A is invertible then W = exp(/ A -1 A') is a matricial 
solution of this system. Hence, the above homogeneous system can be solved at 
precision j by a procedure called HomogeneousResolution in figure^. 

With the same tools, one can check that the following formal expression 
deduced from the formula for variation of constants 

w- 1 J (^(ff) ^ (^.r^-.e.tr)* 

is a solution of system (^|) . This expression can be computed at precision j by 
a procedure called ConstantsVariation in figure |^. 



3.3 Algorithm 

We summarize our algorithm in figure |2|. This is a simplified presentation where 
the technical details are neglected. 

A preprocessing is necessary to construct, from a SLP coding E, another SLP 
which encodes the associated linear variational system VP and the expressions 
used during its integration. This step relies mainly on Theorem 0. 
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The next part of the algorithm consists in the computation at order n + £ + 1 
of the power series solution of VP. We recall that in one iteration, the number 
of correct coefficients is doubled (see Theorem 2 in Q). 

After the main loop, the procedure Coeffs evaluates the SLP VY on the 
series <frj, Tj and Aj where j = hi2(n + £ + 1); this furnishes the coefficients of 



the jacobian matrix (see Section 3.1) 



Last, the rank computations described in Corollary |l| are performed to solve 
the local observability problem. 

If there is more than one output variable, the evaluation of X7Y and the rank 
computations which are necessary to determine (f> can be done in the main loop: 
the computation can be stopped when the expected rank is reached or when 
the computed ranks become stationary. Thus, we can determine the order of 
derivation v and avoid useless computations. 

We now present a rough upper bound for the arithmetic complexity. 
3.4 Arithmetic Complexity Estimation 

Notations: Hereafter, let L denote the complexity of evaluation of the system E 
and let M(j) represent the multiplication complexity of two series at order j + 1. 
Using classical multiplication formula, we have M(j) £ C(j 2 ). 

Furthermore, let N(J) denotes the number of arithmetic operations sufficient 
for the multiplication of two square j x j matrices. Using classical algorithms, 
we have N(n) € 0(f). 

Proposition 2 The number of arithmetic operations on the base field used in 



the algorithm presented in Section 3.S is bounded by 

O \m{v) {N{n + £) + (n + m)L \ +{n + £+ l)N(n + £) 



n + £, 



Proof: From construction done in Section 3.1 and Theorem^, we conclude that 



the complexity of evaluation of the SLP coding dP/d(X,X,Q),'S7P and VY 
is bounded by 0(N{n + £) + (n + m)L). Hence, at each step, the number of 
arithmetic operations necessary to evaluate this SLP on power series truncated 
at order j, is bounded by 0(M (j)(N(n + £) + (n + m)L)) . 

Furthermore, the determination of the first j terms of the solution series of a 
system of linear ODE (||) requires 0(M(j)(N(n) + N(n + £))) arithmetic oper- 
ations by the well-known method of integrating factors (see § 5.2 in ^ for more 
details). So, as M(j) + M{\j/2\) H = 0(M(j)) and as our Newton opera- 
tor is quadratic, the arithmetic complexity of the computations of the jacobian 
matrix 9(£ l G) < i < 1 ,/a(X, 6) is bounded by 0(M(v)(N(n + £) + {n + m)L)). 

To conclude, we notice that the cost of a rank computation for & i x j matrix 
is 0(N(i)j /i) if i < j (see page 108 in Q). The Corollary |l| describes the rank 
computations done at the end of the main loop of our algorithm. ■ 
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Remark: The specialization of input variables on a randomly chosen polyno- 
mial of degree n + £ increases the evaluation complexity L of the system E but 
it does not change the general complexity of the algorithm. When the system 
is not observable we assume that the index v is n + i (see Definition 4 in jTj). 

We have presented the complexity of our algorithm in term of arithmetic 
operations on Q. Such an operation requires a time, roughly, proportional to 
the size of its operands. Using modular techniques, we control the growth of 
the integers involved in the computations. We estimate now an upper bound on 



these integers; this bound will be used in Section 3.6 in order to estimate the 
probability of success of our algorithm. 

3.5 Growth of the Integers 

The forthcoming estimations relies on the formal definition of the jacobian ma- 
trix d(C G)o< i< u I 'd{X , 0) and are not dependent of the computations described 



in Section 3.1 and 3.2 



Let us introduce a measure for the size of a (n + £ 4- r)-variate polynomial 
which influence the growth of the integers (see [g| for more details) . 

Definition 3 Let A be a finite set of non zero integers. The (logarithmic) 
height of A is defined as ht(A) := In \A\ with \A\ := max{|a| + 1, a £ .4}. 

The height of a polynomial with integer coefficients is defined by the height 
of its set of coefficients. 

We summarize in the following lemma some basic properties of height: 

Lemma 1 Let p\, . . . ,p s be (n + £ + r)-variate polynomials with integer coeffi- 
cients, x an integer and d a partial derivation (d/dx for example). 

• ht(dp) < ht(p) + lndegp; 

• ht(p(x j) < ht(x) degp + ht(p); 

• ht(J2i=iPi) ^ max i=i.. « ht(pi) + Ins; 

• ht(pip 2 ) < ht(pi) + ht(p 2 ) + min{degpi,degp 2 } hi(n + £ + r + 1). 



We use the notations introduced in Section 1.1 and we denote by h (resp. d) 
the maximum height (resp. degree) of the numerator and of the denominator 
of the expression involving in system E. 

Proposition 3 Let ho be the maximum of heights of the integers Xq, and of 
the integer coefficients of U . 

• /ii(denom DG (X )) < (2j + l)(n + m) ((2 ln(n + £ + r + 1) + h Q )d + h^j ; 
. hffnumovfirfX W < (^J + ^ + m)((2Mn + £ + r + 1) + h )d + h) 
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Proof: As we are interested in an upper bound, we do not consider the reduced 
form of the fractions fi and gi involved in Cg but we consider that all these frac- 
tions share the same denominator q. So, C — ( fi®i) /<? an d 1 ls the common 
denominator of all gi. 

Thus, the degree of these non-reduced numerators and denominators is 
bounded by (n + m)d and the height by (n + m)(h + d ln(n + I + r + lV) . Let 
us notice that the denominator of C? g is q 2 i +1 ; these facts and Lemma 01 prove 
the first part of our proposition. 

We prove the second part by induction; let us consider (vj)j^n the se- 
quence of polynomials defined by the numerator of g as initial condition vq and 
by the recurrence relation Vj+i :— ^ fi{qdiVj — (2j + \)vjdiq). By construc- 
tion, Vj is equal to the numerator of L? g. Thus, the degree of Vj is bounded 
by (2j + l)(n + m)d — j and we obtain the following recurrence relation from 
Lemma [j]: 

ht(v J+ i) < 2(n + m)(2dln(n + £ + r + l) + h) + ht(vj) +ln2n(2j + l)(n + m)d. 

This is sufficient to conclude. ■ 

Remark: the use of non-reduced fractions simplifies the previous proof but it 
increases the upper bound by a factor (n + m) which is not significant in this 
presentation. 

We have showed that the size of the coefficients of the final specialized ja- 
cobian matrix is mainly linear in the differentiation index v. But some inter- 
mediate computations can require integers of bigger size. In order to construct 
a practical and efficient algorithm, we have to avoid this growth using modular 
techniques. 

Almost all the operations used in our algorithm commute with the canonical 
homomorphism from Q to a finite field F p . But, when we choose a prime num- 
ber p, we have to avoid the cancellation of dPj dX mod t and of the determinant 
of 9(£ l G) < 4 <,/9(X,e). 

The cancellation of dPj dX mod t can be checked at the begining of our al- 
gorithm. Thus, the probabilistic aspects concern mainly the choice of specializa- 
tion and of a prime number such that the determinant of d(C l G)o<i<^ / d(X , O) 
does not vanish modulo p if this matrix is of full generic rank. 



3.6 Probabilistic Aspects 

Hereafter, we call a bad point, a set of specializations {Xq.O,U} where the 
jacobian matrix d(C l G)o<i< v /d(X, 0) is not of full generic rank. Thus, a bad 
point is a zero of the polynomial associated with a minor of this matrix. We 
estimate the probability for a specializations to be a bad point with the following 
proposition. 



Proposition 4 (R. Zippel & J. Schwartz [47]) Let q be a s-variate poly- 
nomial of total degree D and f2 a set of integers. The worst case bound for the 
probability that a point in f2 5 will be a zero of q is D/#Q. 
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This result shows the relation between the choice of the size ho of the used 
specializations and the probability of success of our algorithm. In fact, as 
the determinant of d(£ l G)o<i<v/d(X, 0) is a polynomial of degree bounded 
by D := (n + l){2v + l)(n + m)d, a point in {0, . . . , fi 1 D} ( . n +^ ( . r+1 ~> is not a bad 
point with probability at least 1 — 1/fii- 

Furthermore, we can estimate the probability that the determinant is divisi- 
ble by a prime number p with the following arithmetic analogue of Proposition 0. 



Proposition 5 (§ 18 in [45]) For any integers a and b such that b < a < c, 
the probability that a prime number p between b+1 and 2b divides a is bounded 
by 2 lnc/6. 

From Proposition || and Lemma [j], we can estimate the size of the coefficients 
of the specialization of the jacobian matrix d(£ l G)o<i<„/d(X, 0). Thus, using 
Hadamard's inequality, we find the following rough upper bound for the size of 
the specialized determinant: 

ht(c) := (2 ln(n + £ + r + 1) + h Q )D + (n + t){2v + l)((n + m)h + \n2nD) 

Thus, if the computations are performed modulo a prime number p greater or 
equal to 2ht{c) (12 then the probability that the specialized determinant is not 
divisible by p is at least 1 — l/p2- These results lead to the following estimation. 



Proposition 6 Let n be an arbitrary positive integer and 

D := (n + £)(2u + l)(n + m)d, 
ht(c) := (2\n{n + £ + r + 1) +\nD)D+ (n + i)(2v+ l){(n + m)h + \n2nD) . 

If the matrix d(£ l G)o<i< u /d(X,Q) is of full generic rank then the determi- 
nant of this matrix specialized on random integers in {0, . . . , pD} is not divisible 
by a prime number p > 2ht(c) /i with probability at least (1 — 1//J,) 2 . 



4 Experimental Results 

We present now some benchmarks from an implementation in Maple ]38[ | of 
our algorithm. The Maple computer algebra system provides almost all the 
necessary tools to handle the canonical isomorphism between polynomials and 
polynomial functions: this explains why we have chosen it to implement our 
algorithm. 

The computations summarized in figure || have been performed on a personal 
computer Pentium III (633 Mhz) with 128Mb of memory running Linux 2.2 and 
Maple V.5. This computer was provided by the UMS MEDICIS §§. 

These results show that the index of differentiation is a significant charac- 
teristic of the complexity of algorithm presented in Section |3.3| . Furthermore, 
the last example of the array shows that the complexity of evaluation have a 
significative influence and that the total number of multiplications is clearly less 
significant than the number of multiplications between state and input variables. 
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figure 3: Some benchmarks 



System m v £ n r L time in s. 



V1987 2 8 5 4 17 0.8 

R1986 2 14 9 4 1 19 1.5 

MV1991 2 14 8 5 2 59 2.4 

MW2000 3 18 14 4 67 5.7 

KD1999 2 19 14 5 2 34 6. 

G1995 1 23 17 5 46 10. 

SHH1997 1 23 13 9 38 13.5 



4.1 Certifying the result 

As shown in Corollary [l], the local observability property is associated to the 
fact that the jacobian matrix is of full rank. Our algorithm computes the generic 
rank of this matrix. When it is maximal, the result is certainly correct. Hence, 
if this algorithm states that a model is observable then this result is certified (it 
is a RP-complexity class test, see § 25.8 in |ft5[). 

If there is a non empty set O C X U of non observable variables and pa- 
rameters, the observable parameters can be randomly specialized and there is 
an infinitesimal transformation acting on the non observable state variables and 
parameters, 



which leaves invariant the outputs G and the vector field associated to the 
model. This leads to the following linear system of PDE's: 



This system of PDE can be difficult to solve; nevertheless, we are not interested 
in the whole Lie algebra but in any non trivial subalgebra which can certified 
our result. 

Furthermore, our algorithm decreases the number of unknown of the original 
problem. Hence, in many cases of practical interest, there is a rather straight- 
forward solution (compare with [p0|). For example, these computations have 
been performed in less than a hour with Maple for the following examples. 

4.2 Examples 

We present now the examples indicated in figure [|, the answer of our algorithm 
and some results of the method sketched in section |l.l[ We just give the non 
observable parameters and variables; the other one are observable. 




xeonx 



eeone 




= 0, 
= 0. 
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V1987 Model of a flow reactor to pyrolyze methane 

This example is taken from [Q . 

XI = —X\{kl + k 2 XA) + k 5 X 3 X4, 

±2 = k^xiXi - (k 3 + ki)X2, 

kiX 2 - k 5 X 3 X4, 

Xi(ki + k 2 x 4 ) + 2k 3 X2 — k 5 x 3 X4, 



x 3 

X4 

Vi 
K. V2 



Xi, 
X 2 - 



Our Maple implementation certifies that all the variables and the parameters 
are observable. 



R1986 A pharmacokinetic model 

This example is taken from |32|. The letter u denotes an input. 



XI = U — (Cl + C2)X\, 

±2 = C\Xl - (c 3 + C 6 + Cl)x 2 + C5X4, 

X 3 = C 2 Xl + C 3 X 2 - C4X 3 , 

±4 = C 6 X 2 - C5CC4, 



yi 
m 



CgX2- 



Our Maple implementation gives the following results: 

• the variables {X2, x 3 , X4} and the parameters {c\, C2, c 3 , cj, c%, eg} are not 
observable; 

• the transcendence degree of the field extension k(U, Y) «— > K, is 1. 
Further computations show that the following one parameter group 



x 2 — > Xx 2 

x 3 ((1 - A)c! + c 2 )x 3 /c 2 

X4 — » XX4 

ci — > Aci 

c 2 — » (1 - A)ci + c 2 



c 3 ((1 - A)ci + c 2 )c 3 /Ac 2 

c 7 -> c 7 -c 3 (ci +c 2 )(l - A)/Ac 2 

cs -* -c 8 c 2 /((l - A)ci +c 2 ) 

Cg — > Cg/A 



is composed of symmetries which leave the vector field and the output invariant. 
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MV1991 Model for an induction motor 

This example is taken from j28|. The letters u x and u y denote inputs. 

2/1 = U, 

, 2/2 = * x 2 + *« 2 . 

Our Maple implementation gives the following results: 

• the variables {7 X , I y } and the parameters {M, L s , R s , L r , i? r , J, T{\ are not 
observable; 

• the transcendence degree of the field extension k(U, Y) «— ► K, is 1. 
Further computations show that the following one parameter group 

{I-x,Iyi M, L s , R s , L r , R r , J, Ti} — > {AJx, AJ y , MA, L S A, i? s A, L r / A, Ai? r , AJ, AT;} 
is composed of symmetries which leave the vector field and the output invariant. 

MW2000 Multispecies model for the transmission of pathogens 

This example is taken from J27j . 



b 


= M + ci(j/i +2/12), 








Ai 


= +2/12), 








A 2 


= #2(2/2 + 2/12) + ^2, 








±12 


= (1 - 6»! - 6> 2 )& - (miAi + m 2 A 2 + m)- x i2 + 










(1/1 + r)yi + (u 2 + T)y 2 + ry 12 , 








yi 


= 8xb + miAiZu + v%yi2 - ((1 - 7r 2 )m 2 A2 - 


- f 1 - 


r- M H 


- ci + r)j/i, 


in 


= # 2 fe + m 2 A 2 xi 2 + ^iyi 2 - ((1 - 7Ti)miAi - 


Yv-2,- 


i-ah 


-r)y 2 , 


y\2 


= (1 - 7Ti)miAiJ/ 2 + (1 - 7T2)m,2A 2 yi - (l>i - 


Yv 2 - 


l-M- 


- ci + r)yi2 


0\ 


= %\2 + yi + 2/2 + 2/12, 








0-2 


= 2/1 + 2/12, 








03 


= 2/2 + 2/12- 









Our Maple implementation gives the following results: 
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• with the exception of {f3\, P2, I2, mi, W2}, all the parameters are observ- 
able; 

• the transcendence degree of the field extension k(U,Y) <^-> K. is 2. 
Further computations show that the following two parameters group 

{Pi, fo, h, mi, m 2 } -> {fli/h, P2/I2, h/h, hmi, h^} 
is composed of symmetries which leave the vector field and the output invariant. 

Let us notice that the output o\ is in fact a constraint equal to 1. Hence, 
our model can be composed of relations of order zero which can be considered 
as supplementary outputs. 

KD1999 Model for a chemical reactor 

This example is taken from . 



' Ca 




-C A )-k e- E / RT C A , 


C b 


= -^c B 


+ ko e- E / RT C A , 


t 




-T)-k e- E / RT C A ^ 


To 


= 


T ) U T 3 -T 
3> Ph,c p h V h ' 


Vi 


= Cb, 






= T. 





We denote by A the Arrhenius' law e - E / RT and we add the ordinary differential 
equation A = EAT/(RT 2 ) to the model. Our Maple implementation gives the 
following results: 

• the variable A and the parameters {E, R, AiJ r , U, p, c p , ph, c p h, kg} are not 
observable; 

• the transcendence degree of the field extension k(U, Y) <—* K, is 5. 

Further computations show that the following five parameters group 

A -> XiA 
ko -> fco/Ai 
E — > \%E 
R — > X 2 R 

is composed of symmetries which leave the vector field and the output invariant. 

G1995 Model of Circadian oscillations in the Drosophila period pro- 
tein 

This example is described in introduction. 



P -> X 3 p U _> A3A4C/ 

C P ^4Cp c ph — * XbCph 

AH r -> A 3 A 4 AiJ r Ph -» A 3 A 4j o, 1 /A 5 
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SHH1997 Model of a part of the blood coagulation mechanism 

This example is taken from J40| . 



= kcxX-RVV 
1 kmx+X > 



r 4 



'2 



kpxVa ■ Xa ■ PL, r§ 



ki Xa Xa, 
k PL PT, 



r3 



kc.2,1 1 -Xa 
k'ni2-\-II 5 



f8 



kin aa2 M ■ Ha, r 9 



kcyV -I la 
km v +V ' 

kcn-IIPT 
kmu+II ' 

klllaATIII 



X 


= -ri 








Xa 


= n- 


r 2 - 




+ r 5 , 


V 


= ~T3 








Va 


= r 3 - 


n + 


rs, 




PL 


= -ta 


+ r 5 






PT 


= n - 








II 


= -re 


- r 7 






Ila 


= r 6 + 


ri - 




- r a , 


IIaa 2 M 


= r 9 , 









y 



I la 



556 
1000 



IIaa 2 M. 



Ila. 



Our Maple implementation gives the following results: 

• the parameters {kcx, kmx, key, kmy, kpr, fcc/j, kc 2 } and 
the variables {X, Xa, V, Va, PL, PT} are not observable; 

• the transcendence degree of the field extension k(U, Y) ^ K. is 1. 
Further computations show that the following one parameter group 



X -> 


XX 


PL - 


■+ XPL 


Xa 


XXa 


PT - 


XPT 


V -> 


XV 


kc x - 


-> Xkcx 


Va 


XVa 


kmx 


-> Xkmx 



key 


-» Xkcy 


kmy 


-» Xkmy 


kpr - 


-> k PT /X 2 


ken - 


-> kcn/X 


kc 2 





is composed of symmetries which leave the vector field and the output invariant. 
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